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The phase transition between a massive dense phase and a diluted superparamagnetic phase has 
been studied by means of a direct molecular dynamics simulation. The equilibrium structures of 
the ferrofluid aggregate nucleus are obtained for different values of a temperature and an external 
magnetic held magnitude. An approximate match of experiment and simulation has been shown 
for the ferrofluid phase diagram coordinates “field-temperature”. The provided phase coexistence 
curve has an opposite trend comparing to some of known theoretical results. This contradiction has 
been discussed. For given experimental parameters, it has been concluded that the present results 
describe more precisely the transition from linear chains to a dense globes phase. The theoretical 
concepts which provide the opposite binodal curve dependency trend match other experimental 
conditions: a diluted ferrofluid, a high particle coating rate, a high temperature, and/or a less 
particles coupling constant value. 
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I. INTRODUCTION 

Ferrofluids (FF) are colloidal suspensions of magnetic 
nanoparticles in a carrier liquid. An increasing interest 
to FF is related to their applications: drug deliveries JTj, 
a hyperthermia [2] , a magneto-resonance tomography 0], 
others gH3. Additionally, phase transitions fundamental 
understanding can be developed by means of FF usage 
as a model system. 

Suspended and aggregated phases are two major pos¬ 
sible states of a FF matter [HI E]- A vessel volume con¬ 
tains a single phase or multiple FF phases {20]. A more 
detailed systematization of dense phases is determined 
by an aggregates substructure classification (spatial or¬ 
dering and symmetry) which impacts on thermodynamic 
properties. Some of possible FF aggregates are: drop-like 
aggregates [T0HT2] (either bulk drops or microdrops), one- 
dimensional chains HME2], labyrinthine patterns S3, 
hexagonal patterns [18], rod-shaped [HI l 20] aggregates, 
dumbbell-like aggregates {20] EL], others. Phase tran¬ 
sitions correspond to all possible pairs of these phases. 
Forced phase transitions can be triggered by an external 
magnetic field applying and/or temperature changes or 
through other external impact. 

The phase transition between the suspended and ag¬ 
gregated phases has been called by a liquid-gas (1-g) 
phase transition [91. A contradiction between different 
FF 1-g phase transition research results exists [5]: while 
it is predicted by theories m mm and observed in 
some experiments E3> simulation research reports (e.g. 
!2H!) usually exclude the 1-g phase transition possibility 
or provides a phase diagram which does not match an ex¬ 
pected one [9] . It has been concluded that mean-field and 
statistical models are not justified in the case of large cou- 
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pling constant and/or large densities of FF nanoparticles 
[291 . Other aspects of the above-mentioned contradiction 
will be discussed in the present work. 

A material memory related phenomenon is detected 
experimentally for drop-like aggregates. The phase dia¬ 
gram evolves after a cyclic heating and cooling: temper¬ 
ature and field magnitude transition values depend on 
a sequential cycle number [301 EE]- This effect as well 
as similar ones cannot be explained in scope of theo¬ 
ries which imply a FF as a continuous medium of sep¬ 
arate particles. Such thermodynamical (statistical) the¬ 
ories [22] I25H251I52H55] require assumptions leading to 
an analytical closed-form expression derivation possibil¬ 
ity. Some of these theories are phenomenological (e.g. 
[25] ). Examples of the simplification related assumptions 
of these models are following: the nearest-neighbors ap¬ 
proximation of configurational integral [10] , the isotropic 
potential approximation [35], a dipole-dipole interaction 
as a perturbation {2|] , a very diluted phase consideration 
S3 El, a zero or infinite value of an external magnetic 
field [20], others [39]. Some of these models [36] imply the 
Boltzmann factor based statistical averaging over all pos¬ 
sible magnetic moments pair orientations. However, this 
model assumes an uniform statistical distribution over 
particles mutual positions which is justified only in the 
case of a high-temperature FF with a small correlation 
of magnetic moments orientations [f2Sl . This assumption 
works in the case of a temperature close or higher than a 
critical temperature of a FF. Similar model shows match¬ 
ing an experiment in the case of a diluted enough FF: a 
volume fraction value of a FF dispersed phase was 2.4% 
[24] . This match has been shown only for a sedimentation 
direction of a phase transition in contrast of a “melting” 
[531 of a dense phase. 

A short review and comparison of phase transitions 
investigations in theoretical, experimental and numerical 
simulation research works has been reported in the m- 
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The idealized “homogeneous” FF models treats the phase 
transition as the van der Waals 1-g transition of an ensem¬ 
ble of separate particles. This assumption leads to the 
problems of matching of the theory and the experiment 
m- Any heterogeneous clusters (e.g. one-dimensional 
chains or coil-like structures) are not covered by such 
models. A detailed study of FF aggregates equilibrium 
internal structure shows a much wider variety of possi¬ 
ble structures ITT fo- A thermodynamic stability of FF 
aggregate structures leads to an emergence of a gap be¬ 
tween branches of a binodal of a coexistence of dense and 
dilute phases Ell¬ 
in order to take into account a variety of possible spa¬ 
tial and spin configurations of a dipolar system, a numer¬ 
ical modeling based research work is required. There are 
two major types of such methods: molecular dynamics 
(see e.g. [221 30] and references therein) and Monte Carlo 
(see e.g. [35; ]4T] and references therein) simulations. In 
scope of a conventional assumption that the Brownian 
particles motion does not require a quantum mechanics 
based description [42] , a molecular dynamics is a possible 
candidate to be an ab initio method of a FF numerical 
research. Another assumption is implicit solvation (con¬ 
tinuum solvation) j-131 approximation of Brownian parti¬ 
cle motion. Most molecular dynamics simulations focus 
on systems with periodical boundary conditions HE ED]. 
A comprehensive comparison of a theory, a simulation, 
and an experiment has been reported for magnetization 
curves only [44 1 . Concerning phase diagrams, such com¬ 
parison is either non-quantitative or requires an intro¬ 
duction and a variation of ad hoc unknown parameters 
leading to an expected match: magnetic cores size distri¬ 
bution functions (a polydispersity), a nonspherical parti¬ 
cle shape, the Hamaker constant Ah precise value defini¬ 
tion problem [33], the van der Waals forces microscopic 
theory approximation including many-body interactions 
SSH13, a nonmagnetic part of a nanoparticle volume, 
a stabilizing surfactant layer interaction nature (formula 
and constants) [29] , a coverage rate of a nanoparticle by 
a surfactant, solvation layers [24], aggregating character¬ 
istic parameter [ 251 . etc. There is a lack of numerical 
research works which match experimental phase transi¬ 
tions parameters 1301ED EB3MI even via the selection 
of the above-mentioned ad hoc parameters. According 
to m, existing simulation methods [25] [33 [351 EH 1521 - 
[55] provides mostly linear chain-like aggregates. Hence, 
the phase transition between a massive dense phase m 
and a diluted superparamagnetic phase is not yet enough 
studied by a direct molecular dynamics simulation. 

Our previous published simulation method had been 
designed to the purpose of a monodisperse FF aggre¬ 
gation research in the case of large magnetite particles 
(diameter 20 nm) (T4j. A smaller size of magnetite par¬ 
ticles (~ 10 nm or less) corresponds to a predominance 
of the Brownian motion comparing to aggregating forces 
mm- Hence, in real polydisperse FF complex phase 
transitions and phases coexistence are possible mm- 
The subject of the present research is the phase transi¬ 


tions in the polydisperse FF under applied temperature 
and magnetic field changes. For this purpose, the original 
simulation method should be improved and validated by 
a comparison with related experimental research mm- 


II. MODEL AND SIMULATION 

A. Particle model and interactions 


If a deviation of a nanoparticle shape from a spherical 
one has isotropic probability distribution then one could 
consider particles as hard spheres [56]. The 'i-tli particle 
hydrodynamical diameter d, = d™ + 2ao + 25 consists 
of a ferromagnetic core diameter d™, a nonmagnetic sur¬ 
face layer thickness which equals a lattice constant ao of 
a respective bulk crystalline material, and a stabilizing 
surfactant layer thickness 5 23 E3- We will consider 
only single-domain particles with d™ < do [58] where the 
single-domain particle size threshold is do < 0.5yitm [53- 
m- A stabilizing surfactant molecules surface density Ns 
on a particle determines a rate of its coating and impacts 
on the FF stabilization [56]. 

Lagrangian and Rayleigh dissipation function of a N 
particles system include each particle degrees of freedom 
(i = 1,N): center of mass spatial coordinates r,, Euler 
angles of a particle rotation, and a magnetic moment m; 
where m,; = const. A particle potential energy depends 
on its Euler angles implicitly through an angle 9i between 
m,i and a particle easy magnetization axis unit vector ri; 

m- 

The total force acting on the i-th nanoparticle is given 
by: 


F.-r ;(!) 
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Assuming the limit of low Reynolds number, the vis¬ 
cous friction force is given by the Stokes’ law F- = 
—3nridiVi = — 7 Jvi where the 77 is a carrier liquid dy¬ 
namic viscosity and the 77 is the i-th particle velocity. 
The Fjj d is a dipole-dipole interaction force [H [23 E3- 
The van der Waals’ interaction energy between spherical 
bodies is given in [62] ■ Its differentiation yields the force: 


F w 



-32 A H zRi 3 Rj 3 


(Ri 



— (Rj + Ri) 


( 2 ) 


where Ah is the Hamaker constant; a radius Ri = ao + 
df 1 /2 ; a distance z = (77 — r^j; and its dimensionless 
value pij = [ri — r.j) /z. The energy density of entropic 
repulsion [63] of two surfaces: 


E { = kTN s [1 — s/{2 8)\ 


(3) 








3 


after an integration over surfaces of two particles with different radii takes the form (at z < Ri + Rj + 25): 


Lrij 


irkTNs [z-{2 6 + Rj + R t )f {Rj + R,f {z + 5) - (i ? 3 + ) 


6 6z(Ri + Rj) 


( 4 ) 


where k is the Boltzmann constant; T is a thermodynamic temperature; s is a distance between surfaces; 5 is a 
surfactant molecule length. In the case of particles of same radii, the closed-form expression equation 0 differs from 
the well-known logarithmic one |631 153] . However, these dependencies are numerically close to each other with the 
tolerance ~ 20 % for particles with diameter 10 nm. Existing experimental results cannot justify these expressions 
difference. The comparison of the derivation procedures will be published separately. 

The corresponding force F™ = —dGf^/dr i is given by: 


F% K = -pijnkTNs (z-26-R 2 -R 1 ) {R 2 + Ri) x 


2 (z 2 + 5 2 ) - [R\ - R,R 2 + R 2 ) + l) 


+ 1 + 5z + 1 


2R\R 2 - {R 1 - R 2 f 


R 2 + R\ 
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The hard-sphere model [55] force is given by the following expression: 

pHS f oo • pij , at z 5jj {di + dj )/2 /^\ 

b' ~ y 0, at z > (di + dj)/ 2 ^ ’ 

The total effective torque acting on the i-th particle is determined by a viscous rotational friction, a particle anisotropy, 
the applied external magnetic field B 0 value and a dipole-dipole interaction field 52: '57]: 


n = — 7 » 


K Vi sin (2 di) [rii x m/\ /sin (0,) 


m, x 


Bn 


E b j 


( 7 ) 


where the = 7 rdf 77 ; V t = 7 r(df 1 ) 3 /6 is a i- th particle magnetic core volume; K is a particle magnetic anisotropy 
constant; the u>i is the i-th particle self rotation angular velocity; the Bj is a magnetic field created by the j-th 
particle magnetic moment in the geometric center of i-th particle. All types of magnetic anisotropies should be taken 
into account: a magnetocrystalline anisotropy [65] , a shape anisotropy (demagnetization tensor) [55], an anisotropy 
of surface m , etc. 


B. Translational and rotational motion 


In scope of the continuum solvation approximation of 
not very dense solutions, the Brownian translational and 
rotational particle motion is described by the Langevin 
equations with the hydrodynamic-originated Langevin 
parameters [29l 1451 l52l 1581155] : 


Mid 2 n/dt 2 = F t + g (8) 

IiduJi/dt = Ti + (9) 

where the Mj and /, are the z-th particle mass and mo¬ 
ment of inertia respectively; the t is a time; the and 
are a random force and torque respectively, which are 


usually modelized by Gaussian noise ]2E 40. 00. Sj : 


(€?■(*)> = 0 

(10) 

<X T (t) £.7 (O) = GkftTjJd {t - t') 

(11) 

= 0 

(12) 

(t?(t)Z?(t , ))=6k B T'y?6(t-t') 

(13) 


where <5 (t) is the Dirac delta. 

It is important to note that original Langevin 
equations were supplemented by particles interaction 
forces Q and torques Q. This supplementation has 
been considered as an obvious and intuitive step which 
had been made in scope of the molecular dynamics sim¬ 
ulation based research programs [25, 30, 52j. However, 
it still should be exactly theoretically justified in general 
case of the Brownian particle motion. This statement is 
in need of further research. 

A rotation of a magnetic moment inside a particle is 
described by the Landau-Lifshitz-Gilbert equation with 
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the attempt time to = M s /2a r yK where M s is a satu¬ 
ration magnetization; the a is a damping factor; the 7 
is an electron gyromagnetic ratio m- The simultaneous 
Brownian dynamics of a magnetic moment and a particle 
has been investigated in the m- 


C. Experimental and simulation parameters 

All parameters of the present simulation correspond 
to the specific experiment da eh Eg. The considered 
FF consists of the magnetite nanoparticles suspended in 
the kerosene carrier liquid. The stabilizing surfactant 
is the oleic acid (8 = 2 nm [56]). The particles diame¬ 
ter distribution was determined experimentally [72] . The 
mean diameter is d = 11.5 nm. The magnetite lattice pa¬ 
rameter a 0 ~ 0.8397 nm corresponds to the cubic spinel 
structure with the space group Fd3m (above the Verwey 
temperature) [73], 73] . However, this parameter is ap¬ 
proximate because only 10% of particles have a crystal 
structure. This conclusion was made based on a dark 
field electron microscopy measurements m- The sim¬ 
ulation initial condition is a random close packing m 
(cf. [33]) of particles positions and random magnetic mo¬ 
ments directions. The external magnetic field is uniform. 
Dielectric properties of the surfactant layer are generally 
similar to those of the carrier liquid [55] . Hence, the 
Hamaker constant of a magnetite is Ah = 4 • 10 -20 J 

m- 

A required step of a dense phase emergence is an orig¬ 
inal phase nucleus forming. Consequently, a phase dia¬ 
gram of a nucleus is close to a phase diagram of a bulk FF 
phase. Only the nucleus aggregate will be considered in 
this simulation. Hence, a number of particles N should 
be minimal but not less some threshold where a phase 
diagram start significant changes depending on the N: 
a transition from a bulk phenomenon to a surface one. 
Number N ~ 100 has been selected for the polydisperse 
FF with a lognormal distribution. 

The magnetic moment precession attempt (damping) 
time order of magnitude is tq ~ 10“ 10 — 10 -9 s [56]. In 
the case of particles di < d , the Neel relaxation time 
and the Brownian relaxation time r B relates as [561163| : 

r N < t b < 10 _6 (s) (14) 

Most part of the range 0 < d < d corresponds to the rela¬ 
tion tn << r 0 . In this case both the Neel relaxation flip 
and a magnetic moment dynamics should be considered. 
The Brownian rotation is a much slower process. A parti¬ 
cle rotation does not impact on a magnetic configuration 
and a free energy of the phase. 

The relation is opposite for the larger particles di > d 

m- 

t 0 < 10” 6 (s) < t b < t n (15) 

Here, a magnetic moment alignment with an effective 
field can be modeled as instant. The magnetic moment 


rotates the particle through anisotropy forces. We sup¬ 
pose that the high enough anisotropy energy KV gra¬ 
dient leads to the model with the magnetic moment 
“frozen” into the particle. The vector mi is followed 
by the vector n.^: 

rriiH rii (16) 


The saturated surface density of a number of oleic acid 
molecules in a particle coating is Ag lax = 2-10 18 m -2 [56] . 
Depending on a FF preparation recipe variation (an order 
of mixing/heating of different components |72| . etc.), the 
coating rate k c = N$/N g 10 * can differ. A classical well 
stabilized FF usually has k c ~ 50 % which blocks the par¬ 
ticles aggregation during durable timeframes (years) due 
to the free energy barrier 15 — 25 kT (figure [Taj [SB]. The 
Brownian motion kinetic energy ~ 1 kT is not enough to 
overpass the barrier. Only larger particles form aggre¬ 
gates mm, which corresponds to the potential well at 
the l « 0.5 (figure |Ta|) . In the experimental research [3T] 
the k c value had been reduced, which leads to the dense 
phase (drop-like aggregates) emergence (figure[Tb| . Same 
mesoscopic organization control has been reported in the 
m- The value k c = 5% has been selected for the present 
simulation. The calculation of an entropy and the cor¬ 
responding free energy F by the Eyring’s free volume 
theory m based algorithm [33] has been made for the 
system of particles used in the present simulation (fig¬ 
ure^. The potential well required for the dense phase 
emergence is 2 — 4.5 kT or more 124] . An equilibrium state 
corresponds to the dense phase of such FF. A particles 
contact leads to an infinite negative potential energy of 
the van der Waals interaction. The entropic repulsion 
cannot counteract. However, particles mutual attraction 
is reversible due to an existence of a minimal distance 
between their surfaces s m i n < 8 (figure lb). The value 
Smin = 0.5 ao nnr has been selected as a half of a cell 
constant which qualitatively reflects restrictions of the 
Hamaker theory approximated consideration of spheri¬ 
cal particles as a continuous body ([ 2 ]) in the case of dis¬ 
tances comparable to atomic scale. It corresponds to 
surface structure peculiarities, a nanoparticle quasicrys¬ 
talline structure, etc. On the other hand, the value s m i n 
corresponds to an order of magnitude of two oleic acid 
molecular widths. This additional to the entropic repul¬ 
sion @ stabilized “buffer” role of the surfactant molecule 
has been discussed in the I3B1. 


D. Finite-difference scheme 

The finite-difference Euler scheme is based on the orig¬ 
inal method El- Original method leverages the soft- 
sphere model for a simulation optimization purpose. The 
only changes will be described below. Algorithm details 
of present method can be accessed and contributed in the 
open-source project “Ferrofluid Aggregates Nano Simu¬ 
lator” [73]. The Verlet type of a finite-difference method 
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m is in need of further method improvement. 

A random particle translation and rotation is con¬ 
sidered by means of the viscous limit approximation 
[STHMUHOj which restricts a time tolerance corresponding 
to a finite-difference scheme time step: 

At » max (tJ, t^) (17) 

i v 7 


an inertial term of the motion equations <M can be 
neglected [69] ■ After the integration by the technique 
[US] taking into account ( I0p3 l one can obtain par¬ 
ticle random motion dispersions aj = ^(A?y) 2 ^ and 

af- = ^(A fi) 2 ^ respectively (cf. [29]): 


where tJ = Mi/'yJ and = A/yf". Consequently, 


oj = Dj {6 At + 3 (.Mi / 7 J) [-3 + 4 exp (- 7 ? At / M t ) - exp (-2 7 J At / M z )\ } (18) 

v? = Df{6At +3 (/i/yf) [-3 + 4 exp (-yf At/ h) - exp(-2yf A t / h)]} (19) 


where vector A ipi defines the i-th particle rotation (a 
axis-angle representation; dipi/dt = oij); diffusion co¬ 
efficients are given by the Stokes-Einstein translational 
and rotational equations: Df = kT / 7 ? and Df = 
kT /yb respectively. The corresponding normal distri¬ 
bution (Wiener process of Brownian motion) of the 
and is calculated by the Mersenne twister algorithm 

EH- 


The magnetic moment Neel (14) and Brownian (15) 


relaxation types correspond to different calculation algo¬ 
rithms of the magnetic moment motion which were dis¬ 
tinguished by a particle diameter criterion. In the case 
of the Neel type and At tn, the magnetic moment 
was aligned with the total magnetic field direction in the 
i-th particle center r% at each simulation step. In the 
case of the Brownian type, the rotational motion of 
the particle with the “frozen” (16) magnetic moment is 
calculated. 


III. RESULTS AND DISCUSSION 
A. Equilibration state calculation 

An evolution of the initial random close packing struc¬ 
ture (figure [3]) to the thermodynamically equilibrium 
state was calculated. In the case of the dense phase (the 
nucleus equilibration at the end of the simulation), after 
a stabilization period t s = t s ( Hq,T ), a total moment of 
inertia of the system / to t is not changing except its ran¬ 
dom fluctuations (figure [4]). This statement has been val¬ 
idated up to t = 0.2 s for all obtained dense phases. The 
resulted structures are shown in the figures [5] - [ 6 ] In or¬ 
der to distinguish a diluted and dense phase, simpler and 
similar approach comparing to a conventional correlation 
functions is used. An one-dimensional chain H3H31 is 
defined in the present simulation as a cluster of particles 
which satisfies the condition |/Oj(j+i)| < di + di+i where 
Ath and (i + l)-th particles are neighbors in the chain. 
Here, each particle can have only 1 or 2 neighbors. A 


circle is a particular case of such chain. A set of sep¬ 
arated arbitrary length one-dimensional chains is being 
considered as a diluted phase. If stabilized equilibrium 
state (see figure [4]) corresponds to at least single parti¬ 
cle with the number of neighbors larger than 2 then the 
more complex dense structure is formed. This phase is 
defined as a dense phase by definition. 


B. Phase diagram 

The phase diagram is obtained based on an analysis 
of the resulted structures (figures ©-©> using a bisec¬ 
tion method for the definition of the phase transition 
temperature T t . The binodal curve of the phase dia¬ 
gram (figure [7]) corresponds to the experimental results 
in terms of the trend and the T t approximate value in the 
case of a comparison of temperature values in Kelvins 
[Ml ED- This trend has been observed experimentally 
only starting from the second cycle of the system heat¬ 
ing/cooling. The first heating corresponds to the trend 
weaker than experimental errors and noises 1301 ED- it 
can be explained by the assumption that the initial ex¬ 
perimental sample state does not correspond to the initial 
conditions of the simulation i.e. the aggregate with the 
random close packing m of a polydisperse set of par¬ 
ticles (figure [3]). Indeed, a final equilibrium state of a 
non-stabilizcd FF is a set of primary aggregates which 
consist of large particles only (diameter df 1 ~ 15 — 20 
nm ) M- Hence, the first heating is required in order 
to produce the metastable random close packing struc¬ 
ture which is kept within several heating/cooling cycles. 
The large FF aggregate transforms to the primary ag¬ 
gregates state after durable timeframes which require the 
first heating/cooling cycle again. 

The descending binodal curve (figure [7]) contradicts to 
some theoretical investigations where this dependence is 
ascending [10) [24] (25J. As it is stated at the beginning, a 
contradiction between these theories conclusions and sim¬ 
ulation reports is deeper: even a conceptual possibility of 
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an FF 1-g transition emergence is being discussed among 
different research programs [3]. Mean-field and statisti¬ 
cal models are not justified in the case of a large coupling 
constant and/or large concentration of FF nanoparticles 
[53] . Oppositely, they are justified in the case of a high- 
temperature and/or a diluted phase FF where a corre¬ 
lation in orientations of magnetic moments is small [26] . 
The latter case corresponds to a modeling of a FF by the 
classical Lennard-Jones (LJ) fluid and short-range order 
[35] [551 , which excludes the possibility of an emergence 
of particles complex long-range ordered structures such 
as linear chains, rings, tubes, etc. m ng. However, 
it is well established [3] that an equilibrium FF phase 
microstructure consists of a distribution of chains, rings 
or more complex structures of different lengths. Hence, 
a real FF dense phase should be modeled by a liquid 
crystal like microstructure rather than by the classical 
LJ fluid. However, the LJ fluid is a suitable model for 
conditions which break the long-range ordered structures 
through an entropy increase and a predominance of the 
Brownian motion: a high temperature, a low particles 
concentration (leading to a structures “evaporation”), a 
minor particle coating, etc. In this case a potential en¬ 
ergy summand of a FF free energy [33] is suppressed by 
an entropic summand; the free energy local minimum 
takes place only for large particles (the py < 25 — 30% 
area in the figure [2]) . This is why an 1-g transition in the 
framework of the LJ ferrofluid model is possible only in 
bi-disperse [24] or polydisperse FF thermodynamic the¬ 
ories. 


The LJ like ferrofluid theories correspond to the follow¬ 
ing idealizations: the continuous isotropic phase gaslike 
compression model [25]; a considering of a dipole-dipole 
interaction in scope of the perturbation theory [24] : the 
mean isotropic potential approximation [TUI 1351155] . An 
external magnetic field suppresses the rotational Brown¬ 
ian motion and aligns nanoparticles magnetic moments, 
which leads to a mean attraction force increase [36] . 
Hence, an external magnetic field stimulates the con¬ 
densation phase transition 10] leading to the ascend¬ 
ing binodal curve. In case of a significant correlation 
of magnetic moments orientations (non LJ fluid model 
as discussed above), this conclusion is incorrect because 
the mean isotropic potential model conditions have been 
violated: the significantly anisotropic dipole-dipole inter¬ 
action strongly impacts an equilibrium structure. Such 
structures with a closed magnetic flux (coils, rings, ring 
assemblies, tubes, scrolls, etc. ES!) already correspond 
to a minimal potential energy and a suppressed entropic 
summand. The applying of an external magnetic field 
can only increase the potential energy: flux-closed mag¬ 
netic structures (figure [5a| transform into the set of par¬ 
allel linear chains (figures [5c] 6a 6c I where at least 
marginal particles have a less “workfunction” comparing 
to particles inside the ring or tube structure because they 
are attached to the single magnetic dipole only. Struc¬ 
tures with a closed magnetic flux are bounded by the 
larger dipole-dipole energy. Such aggregate destructs at 


a higher temperature. 

The present simulation results provide a fine structure 
[151115] of the transition from “linear chains” to “dense 
globes” (figures[5]-[6]) through the ring assembly structure 
Hang. Oppositely, the above-mentioned assumptions of 
the theoretical investigations mmm suppress this 
fine structure. The simulation method suggested here 
does not imply idealizations of the theoretical models. 

The last problem which requires discussion is a rela¬ 
tion of the opposite results to experimental studies. Both 
types of binodal dependencies match the respective ex¬ 
perimental observations. Hence, their FF parameters 
were obviously different. The discussed contradiction has 
been shown only for the FF parameters of the present 
research 130 El]: a small particle coating rate k c = 5% 
by an oleic acid which is a concept similar to the [75]. 
This value allows a random compact packing eg in a 
polydisperse FF bypassing a strong surfactant repulsion. 
This is a non-stabilized type of a FF (figure |lb[ ) with a 
high volume fraction of a dispersed phase: an equilibrium 
state corresponds to the value ipy ~ 30%. Oppositely, a 
“classical” stabilized FF (fc c = 50%) corresponds to well 
separated particles (figure |la|). The latter means that a 
correlation in orientations of magnetic moments is small 
[351 due to a larger average distance between particles 
and a corresponding weaker dipole-dipole interaction. In 
this case the aggregates formation can be modeled by an 
1-g type of a phase transition in the framework of the sim¬ 
ple mean field theory [33] or the isotropic potential model 
[351 [353 • In the case of a stabilized FF, an ascending bin¬ 
odal curve theoretical dependence US El Eg matches 
very few experimental reports (the reference in the [24]). 

This research focuses on the FF aggregate nucleus only. 
The phase diagram should be verified and generalized to 
the case of a bulk phase: a larger number of particles 
and a more durable period of a stabilization. According 
to the experiment 130, !3lj, this period should have order 
of magnitude in a range between seconds and minutes 
time-scales. Taking into account an extensive nature of 
the current simulation method m, its natural further 
development is a parallel computing capabilities imple¬ 
mentation. In our opinion, a most promising method 
of the parallel computation of the FF molecular dynam¬ 
ics is a Graphics Processing Units based approach [52] . 
The open-source project dedicated to a further evolution 
of the method has been started [75] ■ A comparison of 
present results with ones based on the Monte Carlo sim¬ 
ulation is in need of further investigations [381, ill] . 


IV. CONCLUSIONS 

1. The equilibrium structures of the ferrofluid aggre¬ 
gate nucleus have been determined for different 
values of temperature and external magnetic field 
magnitude. 

2. The simulation of the ferrofluid phase diagram, 
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which approximately matches the experiment, is 
obtained. 

3. The obtained binodal curve has an opposite (de¬ 
scending) trend comparing to some of known theo¬ 
retical results. 

4. In our opinion, the applied simulation method pro¬ 
vides a more justified description of the magnetite 
ferrofhud with the minor rate of the particle coating 
by surfactant molecules. 

5. The theoretical concepts which provide the ascend- 
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FIG. (1) 

Potential energy dependence on the dimensionless dis¬ 
tance between two identical magnetite particles surfaces 
l = 2 s/d. Different particles diameters d are specified. 



<Pv [ % ] 

FIG. (2) 

Average free energy of the particle in the aggregate with 
100 particles with random positions which fits to the den¬ 
sity. The dependence on the volume fraction of the fer- 
rofluid dispersed phase is shown. Coating rate by the 
oleic acid is k c = 5%. Different particles d diameters are 
specified. 



FIG. (3) 

Initial conditions: perspective projection. A space is lim¬ 
ited by the cubic vessel with the edge length L = 0.3 ^m. 
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FIG. (4) 

Total moment of inertia time dependence at the condi¬ 
tions: the temperature T = 25 °C and the external field 
Hq = 0 Oe. The obtained approximated value of the 
dense phase stabilization period is t s ~ 5 ms. 



a: Ho = 0 Oe and T = 
50 °C 



c: Ho = 200 Oe and T = 
35 °C 




d: H 0 = 200 Oe and T = 
55 °C 


FIG. (5) 

Resulted thermodynamically equilibrium structures: 
perspective projections (area with largest particles is 
highlighted); L = 0.3 /iin. A size of volumetric arrows is 
proportional to a geometrical size of the corresponding 
nanoparticle for an illustrative purpose only. 
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a: H 0 = 400 Oe and T = 
15 °C 



0 °C 


b: H 0 = 400 Oe and T = 
20 °C 



FIG. (6) 

Resulted thermodynamically equilibrium struc- 
tures(continue). 



FIG. (7) 

Ferrofluid nucleus phase diagram: a region “A” - the ag¬ 
gregated phase and the diluted phase coexistence in the 
vessel volume; a region “B” - the diluted superparamag- 
netic phase only. The simulation results (circles) form a 
binodal curve. The experimental results l-'S 1 [ (diamonds) 
are given for a comparison. 





















